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ABSTRACT 

Optical and infrared emission lines from Hll regions are an important diagnostic used to study 
galaxies, but interpretation of these lines requires significant modeling of both the internal structure 
and dynamical evolution of the emitting regions. Most of the models in common use today assume 
that H II region dynamics are dominated by the expansion of stellar wind bubbles, and have neglected 
the contribution of radiation pressure to the dynamics, and in some cases also to the internal structure. 
However, recent observations of nearby galaxies suggest that neither assumption is justified, motivating 
us to revisit the question of how H II region line emission depends on the physics of winds and radiation 
pressure. In a companion paper we construct models of single Hll regions including and excluding 
radiation pressure and winds, and in this paper we describe a population synthesis code that uses 
these models to simulate galactic collections of Hll regions with varying physical parameters. We 
show that the choice of physical parameters has significant effects on galactic emission line ratios, and 
that in some cases the line ratios can exceed previously claimed theoretical limits. Our results suggest 
that the recently-reported offset in line ratio values between high-redshift star-forming galaxies and 
those in the local universe may be partially explained by the presence of large numbers of radiation 
pressured-dominated Hll regions within them. 

Subject headings: galaxies: high-redshift — galaxies: ISM — HII regions — ISM: bubbles — ISM: 
lines and bands 



1. INTRODUCTION 

Ratios of optical and infrared lines from Hll regions 
are popular diagnostics that have been used to infer 
a large number of properties of galaxi es. Perhaps the 
most famous example of this is the 'Baldw in et al. 
( 1981| ) diagram (hereafter the BPT diagram), which 
pI5tr[0 iii]A5007/H/3 versus [N ii]A6584/Ha. Hll re- 
gions in the local universe form a narrow sequence in 
this diagram, and their position along this sequence pro- 
vides information about properties of the Hll region 
such as its density and metallicity. Rec ently, thanks 
to the Slo an Digital Sky Survey (SDSS) ( [Brinchmann 
|et al. 2QQ4| [Tremonti et al. |2004 ) the sequence has been 
extended to unresolved galaxies in the local universe. 
The SDSS showed that galaxies whose line emission is 
dominated by an active galactic nucleus (AGN) or by 
fast shocks are distinguishable in the BPT diagram from 
those whose emission is powered predominantly by star 
formation. Star-forming galaxies and H II regions in the 
local universe follow the same sequence, suggesting that 
star forming galaxies can be simplified as a collection of 
Hll regions. In contrast, AGN-dominated galaxies lie off 
this sequence. 

However, star-forming galaxies at high redshift appear 
to be offset (upward and to the right) i n the BPT di- 
agram fro m thos e in the local u niverse ( Shapley et al. 
2QQ5[ pu et al. 20081 |Erb et al.|[2Q06^ ,2010] ), but do not 
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occupy the same locus as local AGN-dominated galaxies 
either. Several possible causes for the offset have been 
suggested. One possibility is that Hll regions at z ~ 2 
follow the same star-forming sequence as in the local uni- 
verse, but the presence of an unresolved AGN or shocked 
gas contaminates their l ine emission cau sing the shift in 
the BPT diagram (e.g. 



Liu et al. 



ghi 



2008D. Observational 



support for this idea comes from Wright et al. (2010), 
who demonstrate using integral field spectroscopy that 
a weak AGN is resp onsible for the shift of a z = 1.6 
galaxy. Trump et al.| (|2011) stack HST grism data from 
many galaxies to show that this phenomenon is reason- 
ably common. However, another possible explanation 
for the offset is that there are systematic differences ex- 
ists between H II regions in the local universe and at high 
redshift. This suggests that the time is ripe for a reinves- 
tigation of the physics driving Hll region line emission, 
and thus the location of galaxies in diagnostic line ratio 
diagrams such as the BPT plot. 

The problem of computing the integrated line emission 
produced by a galaxy containing many H II regions can 
be roughly decomposed into two separate steps. The first 
is determining the internal structure of an Hll region 
given its large-scale properties, for example the radius of 
the ionization front and the luminosity of the star cluster 
that powers it. The second is determining the dynamics 
of the H II region population in a galaxy, which sets the 
distribution of H II region properties. The first of these 
problems is generally solved with by a radiative transfer 
and chemic al equilibrium co de such as Cloudy ( Ferlandl 
let a m998|) or MAPPINGS ([Sutherl a nd fc Dopita ,^ 1993j 
Dopita eral.|2000[ |Kewley et al.||200l| ), while the second 
is solved by a population synthesis code that generates a 
population of H II regions and follows their expansion in 



2 



the interstellar medium (e.g. Dopita et al. 2006b). For 
this second step, the results depend on what drives Hll 
region expansion, i.e. whether Hll regions are classical 
Stromgren sp heres whose expansion is driven by warm 
gas pressure ( Spitzer| |1978), wind bubbles whose expan- 
sion is controlled by tn e pressu re of shocked s tellar wind 
gas (Castor et al. |l975 ; Wea ver et al.|p!977t, radiation 



pressure- driven sneTTs (Kr umholz &: Matzner|2009 Mur 
ray et al.]|2010 ), or something else. 



I'he m ost commonly - used popula tion synthesis models, 
those of jgqpita et"aL| ([20001 |2005l |2006a b ) and Groves 
et al. ( 2008p , assume that the expansion of Hll regions 
IS prmiarily wind-driven. However, recent resolved ob- 
servations of H II regions in nearby galaxies ha ve shown 
that this assump tion is likel y to be incorrect. IHarper-| 
Clark fc Murray! (|2009| and (L opez et a l.| (|201lj use X- 
ray observations of Carina and 30 Doradus, respectively, 
to directly estimate the pressure of the shocked hot gas 
inside expanding H II regions]^ By comparing these pres- 
sures to the other sources of pressure driving the expan- 
sion, and to the values expected for a wind bubble so- 
lution, they conclude that the giant Hll regions cannot 
be expanding primarily due to shocked wind gas pres- 
sure, and that radiation pr e ssure may well be dominant. 



Moreover, Yeh & Matzner (2012) found no evidence for 



wind-dominated bubbles either in individual regions or 
on galactic scales, using observed ionization parameters. 
Physically, the surprisingly weak role of winds is likely a 
result of Hll regions being "leaky", so that the hot gas 
either physically escapes, or it mixes with cooler gas, and 
this mixi ng; cools it enough for radiative losses to become 
efficient (Townsley et al.|2003). Regardless of the under- 



lying cause, though, the observations clearly show that 
the wind bubble model should be reconsidered. 

In this work, we investigate the implications of these 
observations, and more broadly of varying the physics 
governing Hll regions expansion, for line emission and 
line ratio diagnostics. To do so we create a popula- 
tion synthesis model in the spirit of Dopita's work, and 
within this model we systematically add and remove the 
effects of radiation pressure, and w e vary the stellar wind 
strength. In a companion paper (|Yeh et al. 20121 here- 



after Paper I) we generate a series ot hydrostatic equi- 



librium models o f Hll region s using Starburst99 ( Lei 
therer et al.| 1999|) and Cloudy ( [Ferland et al.|1998[ ), botli 



including and excluding radiation pressure and stellar 
winds. In this paper we use these models to predict the 
integrated line emission of galaxies containing many H II 
regions. 

The structure of the remainder of this paper is as fol- 
lows. In § [2] we describe the method we implement to 
generate synthetic galaxies. In §[3] we analyze the main 
results, with particular attention to how various physi- 
cal mechanisms affect observed line ratios, and in §[4] we 
compare to observations. We finish with discussion and 



^ Note th at|Pellegrini et al.| ( |201l| analyze the same region (30 
Doradus) as'Lopez et al.l(l'2UIl| and report a much higher pressure 
in the X-ray emittmg gas, such that this pressure exceeds radiation 
pressure. They reach this result by ado pting a small filling factor 
for the X-ray emitting gas, compared to [Lopez et al.[ s assumption 
of a filling factor close to unity. However, with such a small filling 
factor, the hot gas is not dynamically important for the H II region 
as a whole, and thus the general conclu sion that hot gas is dynam- 
ically unimportant remains true even if |Pellegrini et alT| s preferred 
filling factor is correct. 



TABLE 1 
Static H II region Models 



Model 


Prad 


\ogn 


RPWW (Radiation Pressure Weak Winds) 


yes 


-1.5 


RPSW (Radiation Pressure Strong Winds) 


yes 


2 


GPWW (Gas Pressure Weak Winds) 


no 


-1.5 


GPSW (Gas Pressure Strong Winds) 


no 


2 



conclusions in § [5| 



2. METHOD 



We are interested in the computing the total line emis- 
sion of multiple H II regions in an unresolved galaxy, such 
those at high redshift, in order to create a synthetic set 
of data that is directly comparable with observed galax- 
ies in the BPT diagram or similar line ratio diagrams. 
The procedure consists of two parts. First, we create 
synthetic line emission predictions for a variety of sin- 
gle Hll regions over a large grid in stellar luminosity, 
radius, and age. We describe this procedure in detail in 
Paper I, but for convenience we briefly summarize it be- 
low. Second, we build a population synthesis code that 
creates, evolves and destroys Hll regions, and computes 
the summed line emission. 

2.1. Spectral synthesis and photoionization models 

We create a population of static, single Hll regions, 
with a wide range of sizes and ionizing lumino sities. To 
do so we use Starburst99 ( Leitherer et al.|1999 ) to gener- 
ate ionizing continua from coeval star clusters of different 
ages. We feed the syntheti c spe c tra in to Cloudy 08.00, 
last described by Ferland et al. (|1998l) , as the ionizing 
continuum emitted at the center oTeach simulated H II 
region. Each H II region is spherical and in perfect force 
balance. We adopt Cloudy's default solar abundances 
and ISM dust grain s ize distributio ns, and the same gas 

We compute 



phase abundances as Dopita et al. 
a grid of models covering a wide range in density, from 
log ^H, in = —1 to 5, where nH,in is the number density 
of hydrogen nuclei at the inner boundary of each H II re- 
gion. Each set of the simulations outputs the integrated 
luminosity of selected optical emission lines, including 
Ha, H/3, [O lll]A5007, and [N ll]A6584, i.e. the lines that 
enter the BPT diagram. For more details we refer readers 
to Paper I. 

We compute four sets of static H II region models, cor- 
responding to four combinations of radiation pressure 
(Prad) and stellar wind strength (Table [T]). In the models 
with Prad 7 radiation pressure is allowed to exceed ionized 
gas pressure, in contrast to Cloudy's default setting. For 
models where radiation pressure is absent, the outward 
force due to the incident radiation field is turned off. We 
parameterize the strength of the stellar wind by which 
is defined as 

P V 

n = ^ S'"" , (1) 



PlF^IF 



Pin Viri 



where Pif'^if — Pin^in is the difference of the product 
of gas pressure and volume between the ionization front 
(Pif'^f) and the inner edge of the Hll region (PinMn), 
which is the outer edge of a hot, wind-pre ssurized bubble. 
Q is th e same wind parameter defined in | Yeh fc Matzner] 
(2012), and we refer readers to Table 1 and Section 4.1 in 



3 



TABLE 2 
Fiducial parameters 



Parameter 



Value 



V[a,,min 


20 M0 


Ma,,max 


5x10^ Mq 


kp 


1 


M 


30 




10^ K cm-3 


SFR 


1 M0 yr-i 


Prad 


yes 


/trap 


2 




0.73 




3.2 



that paper for detailed discussion of its meaning. How- 
ever, an intuitive explanation of Vt is that it measures 
the relative energy content of the hot stellar wind gas 
and the warm photoionized gas; high values of corre- 
spond to wind-dominated Hll regions, while low values 
to ones where winds are dynamically unimportant. 

2.2. Population synthesis code 

We treat a galaxy as a collection of H II regions only, 
with no contribution to line emission from other sources 
(e.g. stars or warm ionized medium). We generate, 
evolve, and destroy these H II regions using a population 
sy nthesis code derived fro m t he GMCEVOL code descr ibed 
in |Krumholz et al.| ( |2QQ6l ) and |Goldbaum et al.| ( [2QTT] ). In 
our models, we characterize a galaxy by two parameters: 
a (constant) star formation rate (SFR) and a mean am- 
bient pressure Pamb, and we give fiducial values of these 
parameters in Table [2j though below we explore how our 
results depend on these choices. For all the results de- 
scribed in this paper, we run our simulation code for 200 
Myr, and write output every 1 Myr. We describe each 
step the code takes below. 

Creation — To create Hll regions, we pick a series of 
stellar association masses Ma from a probability distri- 
bution 

p{Ma) (X (2) 



in the range Ma,min to Ma,max (Williams & McKee 1997). 
We give fiducial values of the mmimum and maximum 
masses in Table [2j but experimentation shows that these 
choices have almost no effect on our final result. Each 
association appears at a time dictated by the SFR; for 
example, if the first three associations drawn in a cal- 
culation have masses of 10^ M©, 10^ M©, and lO'^ M©, 
and the SFR is 1 Mq yr~^, the first association turns on 
at 0.1 Myr into the simulation, the second at 1.1 Myr, 
and the third at 11.1 My r. When an association turns 
on, we pick stars from a |Kroupa| ([2001) IMF until we 
have enough stellar mass to add up to the association 
mass. For computational convenience we discard stars 
with masses below 5 M©, since these contribute negligi- 
bly to the ionizin g luminosity. For the stars we retain, 
we use the fits of Parravano et al.| ([2003) to assign an 
ionizing luminosity and a main sequence lifetime. Each 
association becomes the power source for a new H II re- 
gion, with an ionizing luminosity determined by the sum 
of the ionizing luminosities of the constituent stars. Note 
that we account for aging of the stellar population in the 
ionizing spectrum, but use step-function approximations 
for the luminosity, ionizing luminosity and wind trapping 



factor in our dynamical calculations. These we take to 
be constant during the ionizing lifetime of each cluster. 

Expansion — The neutral gas in which each Hll region 
expands has a radial density profile p = Po(r/ro)~^p, and 
our code allows /c^, = or 1. As we discuss in Section|3.3 



this choice proves to make very little difference, so unless 
stated otherwise we simply adopt kp = 1. We determine 
the mean values of po and ro from two constraints, one 
related to the pressure of the galaxy and a second from 
the mass of the association. Specifically, we require that 



Ma = [47r/{3-kp)]pofl 



(3) 
(4) 



The first of these equations is equivalent to the state- 
ment that the mass of the association is comparable to 
the mass of the surrounding gas (i.e. that the star forma- 
tion efficiency in the vicinity of an association is ~ 50%), 
while the second is equivalent to the statement that the 
gas around an association is in approximate pressure bal- 
ance with the mean pressure of the galaxy. These two 
statements uniquely determine pofo, but we add a ran- 
dom scatter on top of this to represent the expected den- 
sity variation present in a turbulent medium. Such media 
have density distributi ons well- described by lognormal 
distributions (e.g. Padoan fc Nordlund^^2002 ). We there- 
fore scale our value of po by a factor x drawn from the 
distribution 



p{x) 



1 



v/2^ 



exp 



(in X — \nx) 
2^ 



(5) 



where \nx = the dispersion of pressures is Gx = 

Y^ln(l + A1^/4), and M is the Mach number that char- 
acterizes the turbulence. Thus the final value of po^o 
we adopt for a given H II region is p^f^x, with pn^o de- 
termined by the solution to equations ([3| and (p|, and 
X chosen from the distribution given by equation ([5|. 
We adopt a fiducial Mach number = 30, appropri- 
ate for giant molecular clouds in nearby galaxies, but we 
have experimented with values up to = 3 00, appro 



priate for ultra lum inous infrared galaxies (see Krumholz 
Thompson|2007| for more detailed discussion) . We find 



that the choice oi M. makes little difference to the final 
result. 

Once we have the density distribution around an H II 
region, we can compute its expansion. We do so in two 
possible ways . The first is simply following the classical 



Spitzer ( 1978 ) similarity solution for gas pressure-driven 
expansion, generalized to our density profile. The second 
is using the Krumholz & Matzner (2009) generalization 
of this solution to the case where radiation pressure is 
dynamically significant. For this case, we u se the ap- 



Krumholz 



proximate s olution given by equation (13) of ^ 

|fc Matzner l This solution involves a few free parame 
ters, and the values we adopt are summarized in Table 
[2j The most important of these is /trap, which repre- 
sents the factor by which trapping of photons and wind 
energy within the expanding dust shell amplifies the ra- 
diation pressure force. We adopt a relatively low value 
/trap = 2 as a fiducial value, based in part on recent 
simulations indicating the r adiative trapping is likely to 
be very inefficient (Krumholz fc Thompson^2012a|b| , but 
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we also explore different v alues of Aran below. No t e that 
in the case /trap = 0, the Krumholz fc Matzner (2009) 
solution reduces to the classical ^:^pitzer Xl978| ) one. We 
discuss the remaining free parameters below. Finally, 
note that we do not consider the case of expansion fol- 
lowing a Weaver et al. (1977) wind bubble solution, both 
because Dopita et al.| ( [2006b ) have already obtained re- 
sults in this case, and because the observations discussed 
in the Introduction suggest that this model is unlikely to 
be correct. 

Stalling — We stop the expansion of an Hll region if its 
internal pressure ever falls to the pressure of the ambient 
medium (Pamb)- We can express the internal pressure 
as the sum of the thermal pressure of the ionized gas 
and the radiation pressure. The thermal pressure of the 
ionized gas Pgas is 



Pgas = MnniimHCn 



(6) 



where nu is the number density of hydrogen nuclei in 
the Hll region, cn = 9.74 km s~^ is the sound speed, 
/iH = ^/X = 1.33 is the mean mass per H nucleus in units 
of amu, and X = 0.75 is the hydrogen mass fraction. We 
derive nu from photoionization balance, which requires 
that 



(l)S = -Trr^aBniiTie 
o 



-Trr^aB 



Y 
4X 



(7) 



where S is the number of ionizing photons per second 
injected into the region, nu is the number density of H 
nuclei, Ue is the number density of electrons and is the 
case-B recombination coefficient. The factor l-^-Y/AX = 
1.1 (assuming helium mass fraction Y = 0.23, and that 
He is singly ionized) accounts for the fact that there are 
electrons from He as well as from H, and the factor of 
(j) = 0.73 accounts for ionizing photons that are absorbed 
by dust instead of hydrogen. Thus we have 



nil 



3(I)S 



(8) 



Note that this expression implicitly assumes that the 
density within the Hll region is constant, which is not 
the case if radiation pressure exceeds gas pressure. How- 
ever, in this case the gas pressure is non-dominant, so it 
matters little if we make an error in computing it. The 
radiation pressure Prad is 



rad 



trap 



47rr^c 



(9) 



where ip = L/{Seo) is the ratio of the star's bolometric 
power to its ionizing power counting only an energy eo = 
13.6 eV per ionizing photo n . We adopt ^0 = 3 .2 fol lowing 
Murray fc R ahman (20101), iFall et al.l (|2010|) , and iLopez 
let al.lpOllI). 



Destruction — We remove an Hll region from our calcu- 
lation when the stars that provide half its total ionizing 
luminosity reach the end of their main sequence lifetimes. 
This may occur before or after stalling, depending on the 
ambient conditions. 




1 10 100 
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Fig. 1. — Radius of the ionization front (Rw) versus ionizing lu- 
minosity (S). Each data point represents a single HII region at 
one time step of one of our population synthesis calculations; the 
line against which points have accum ulate d on the right side of the 
plot is the stalling line (see Section \2.2\ . Colored lines indicate 
radii and ionizing lumi nosit ies of the M il regions in the RPWW 
model grid (see Section [2J^ ; note that only a subset of the models 
are shown in order to avoid clutter. Blue colors indicate mod- 
els with constant density log nn, in 5 and red colors indicate models 
of constant logi^ip/ ^ch, where r^h is the characteristic radius at 
which radiation and gas pressure balance [Yeh &l Matzner] ( |2012| ); 
note that fch is a function of S alone, and does not depend on 
RiF- The values of lognH,in and logi^ip/ ^ch are as indicated in 
the color bars. 

2.3. Calculation of the line emission 

The population synthesis code generates output files 
containing information about the H II regions present at 
each timestep. For each Hll region, we keep track of the 
ionizing luminosity S of the driving stellar association, 
the radius Rip of the ionization front, and the age t of the 
association. In order to assign line emission luminosities 
to each H II region, we perform a three-dimensional inter- 
polation on i?iF, and t, using the tables of individual 
H II region models described in Section |2.1[ 

Figure [l] illustrates the procedure. The Figure shows 
the ionization front radii (Rw) and ionizing luminosities 
{S) of all the Hll regions present at a single time step 
in one of our population synthesis calculations, overlaid 
with a grid of models for single H II regions at an age of 
Myr. The model grid is characterized by values of density 
^H,in at the inner edge of the Hll region and by the ratio 
of the ionizatio n front radius to the cha racteristic radius 
fch, defined by Yeh fc Matzner' ('2012') as the value of 



RiF for which gas pressure and unattenuated radiation 
pressure at the ionization front are equal. This radius is 
given by 

~ ^^^^ nn^ 

12^(2.2Mic)2^' ^ ' 

where /c^ is the Boltzmann constant, L is the bolometric 
luminosity, Tn = 7000 K is the temperature of the ion- 
ized gas and the factor 2.2 is obtained by assuming that 
helium is singly ionized everywhere. Since L = V^5'eo, the 
value of fch is simply proportional to S. For the simplest 
case of H II regions with an age of Myr, we assign each 
one a luminosity in the [O III] , [N ll] , Ha and H/3 lines by 
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single Hllreg • 
integrated galaxy + 



1.0 -0.5 
log([NII]/Ha) 



0.0 



Fig. 2. — Result of the interpolation plotted on the BPT plot for 
one snapshot of a simulated galaxy. The dashed line is the'Kauff-j 
[mann et al. (2003) line, an empirical separation of star forming 
galaxies and AGIN based on the SDSS galaxies. The dotted line is 
the pKewley et al. (2001) theoretical limit for star forming galaxies. 
Thedot-dashed line is a fit to the star forming galaxies from the 
SDSS galaxies Brinchma nn et^ al.| ( |2008| ). Each HII region is plot- 
ted with a dot (as in l:^ igure[T]) and tiie integrated galaxy is shown 
with the plus sign. 

interpolating between the line luminosities of the near- 
est points in the overlaid model grid. The procedure for 
older Hll regions is analogous, except that there is an 
additional interpolation in age. Once we have assigned 
a luminosity to each Hll region, the total line luminos- 
ity of the galaxy is simply the sum over individual Hll 
regions. Figure |2] shows the final result, where we have 
used the computed line ratios of both the individual H II 
regions from Figure [l] and the integrated galaxy to place 
them in the BPT diagram. 

3. RESULTS 

Our aim is to investigate how radiation pressure and 
stellar winds affect galaxies' emission line ratios. As dis- 
cussed above, the effects are both internal - changing 
the density distribution and thus the emission produced 
within single H II regions - and external - changing the 
distribution of H II region radii and other properties. It is 
easiest to understand the results if we tackle t he in ternal 
effects separately first, which we do in Section 3.1 Then 
in Section [3^ we consider external effec ts an d how these 
interact with internal ones. In Section 13.31 we consider 
how the results depend on the properties of the galaxy 
as a whole (e.g. star formation rate, ambient pressure). 

3.1. Internal effects of radiation pressure and winds 

We first examine how our four internal structure mod- 
els from Table [l] distribute Hll regions in the BPT dia- 
gram. 

3.1.1. Models with weak winds 

We compare the two models with weak winds, RPWW 
and GPWW, in Figure |3] We show H II region models 
with constant lognH,in = — 1 to 5 (blue) and with con- 
stant log RiF I Teh (red) , where R\y is the ioniz ation front 
radius and fch is the characteristic radius in |Krumholz| 





1.5 




1.0 




0.5 




0.0 




-0.5 




1.5 




1.0 


X 


O 


0.5 


D) 




O 


0.0 




-0.5 




1.5 




1.0 




0.5 




0.0 




-0.5 






\ 


\\ \ ■ 

/. ... . 1 .. 


2 Myr 




2 Myr 

\ \ 
\ \ 






\ 

^ \ '\ 

; \ 



-1.5 -1.0 -0.5 0.0 -1.5 -1.0 -0.5 0.0 

log([NII]/Ha) 

Fig. 3. — BPT diagram for the models with weak winds evolving 
from to 4 Myr (from top to bottom). The left column shows 
Model RPWW (with radiation pressure) and the right column 
shows Model GPWW (without radiation pressure). The models 
are shown with lines of constant lognH,in = — 1 to 5 (blue) and 
constant logi^ip/fch (red), where K\y is the ionization front ra- 
dius and rch is the characteristic radius in |Yeh Matzner| ( |2012| > 
at which radiation and gas pressure balance. 



'fc Ma tzner] ([2009'') at which radiation and gas pressure 
balance. The ratio R\y / Teh is related to the ionization 
parameter, as discussed in Paper 1. Within each column, 
we plot three stages of the evolution of the cluster: 0, 2 
and 4 Myr (from top to bottom). 

We confirm some trends that have been seen in the past 
dPopita et a"L]l2QQQ[ |2QQ6b[ [Kewley et aT||2QQlD , such as 
the decrease of line ratios as the cluster ages and increase 
of the ionization parameter from bottom right to top left. 
We explore for the first time a large range of values for 
the density. We find that the higher the density the 
stronger the [N ll] and [O ill] emission, up to the point 
that the gas density exceeds ~ 10^ cm~^. Beyond this, 
the density in the HlI region exceeds the critical densi- 
ties of the [N ll] and [O ill] lines (6.6 x 10^ cm and 
6.8 X 10^ cm~^, respectively) causing the line intensity 
to stop increasing. However, before this point is reached, 
in the highest density models the [O ill] emission is large 
enough that the [O lll]/H/3 ratio exceed s the upper limit 
for starburst models described by [Kewley et al. (2001) 
(black dotted line in Figure [3| . 
We can unde rstand why our n iodels exceed the |Kew- ' 



ley et al.| (2001 ) limits as follows. Kewley et al. I created a 
grid of photoionization models with tixed initial density 
nH,in = 350 cm~^ and strong stellar winds (i.e. assum- 
ing planar geometry), with a range of metallicities and 
ionization parameters, and without the effect of radia- 
tion pressure. They find that the line ratios in their 
model never exceed the limit indicated by the black dot- 
ted line in Figure [3) Our models exceed this limit be- 
cause they reach regimes of very high density and very 
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high radiation field that the |Kewley et al.| models, due 
to their assumption of a fixed density and planar ge- 
ometry, are unable to access. The underlying physical 
processes become clear if we compare our various mod- 
els. Both mo dels GPWW and RPWW can exceed the 
[Kewley et al. limit, while our strong wind mo dels ei ther 
do not exceed or barely exceed it (see Section 3.1.2). In 



model GPWW, the density of the gas near the ionizing 
source can remain unphysical high even when the lumi- 
nosity is very high; as a result there is significant emis- 
sion from high-density, highly-irradiated gas. By con- 
trast, in model RPWW, strong radiation pressure pushes 
gas away from the ionizing source when the luminosity 
is high, which in turn reduces the amount of gas that is 
both dense and highly irradiated. This model still breaks 
the [Kewley et al. limit, but by less than GPWW. When 
stellar winds are included, on the other hand, the wind 
pushes the gas away from the source, reducing the radia- 
tion flux it experiences. This strongly limits the amount 
of dense, highly-irr adiated gas in b oth of our strong wind 
models, and in the [Kewley et aL] models. We therefore 
see that the Kewley et al. limit is not a limit imposed 
by the physic ~of Hll r egions in general; instead, it is 
driven by Kewle y et al.[ s assumptions about the struc- 
ture of Hll regions, and the limitations on density and 
ionizing luminosity that these assumptions imply. 

Comparing the cases with and without radiation pres- 
sure, we see that models with radiation pressure often 
produce less [O ill] emission that those without. This 
effect arises because H II regions with radiation pressure 
and that have Rw / ^ch ^ 1 have most of their gas in a 
radiation-confined shell that has a steep density gradient. 
This should be compared to the mostly uniform density 
produced if one i gnores radiation pressure (Draine 2011; 



Yeh & Matzner 2012; Paper I). The higher density in 
this shell means that the density in the bulk of the emit- 
ting gas can exceed the critical density for a line even 
when the mean density of the Hll region is below this 
value. Hence, Model RPWW saturates at a lower value 
of [O 111] /Up than model GPWW. 

3.1.2. Models with strong winds 

We show the BPT diagram locations of the strong 
stellar wind models, RPSW and GPSW, in Figure [4j 
The first thing that is evident from the Figure is that 
Model RPSW produces line ratios in the BPT diagram 
far from both the other models and from the locations 
of observed galaxies. The region of parameter space 
where the models are not physical within the context 
of RPSW corresponds to Hll regions with large ion- 
izing luminosities but small radii, and one can under- 
stand why Model RPSW avoids this region with a small 
thought experiment. A value of Q = 100 implies that 
PifVif / PinVin — 1 <C 1 (scc Equatiou IT]) , meaning that 
the shocked wind gas dominates the total energy budget. 
This in turn requires that Vw ^ Vin and Pip ~ Pin, so 
that the wind bubble fills almost the entire volume of the 
H II region, leaving only a thin shell of photoionized gas, 
and the gas pressures are nearly identical at the inner 
and outer edges of this shell. However, if the ionizing lu- 
minosity is large enough (specifically if it is large enough 
so that RiF < ^ch), this is impossible. As 5 ^ oo the 
radiation pressure at the inner edge of the photoionized 
shell must greatly exceed the gas pressure, and the gas 
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Fig. 4. — Same as Figure [3] but for the strong-wind Models 
GPSW and RPSW. 

pressure Pjp at the outer edge of the H II region, where 
all of the radiation has been absorbed, must be equal to 
the total pressure at the inner edge, which is the sum of 
the small gas pressure Pin and the much larger radiation 
pressure. It therefore follows that at sufficiently large S 
one mus t have Ptf / Prad ^ 1 , ^ ivin^ 17 <C 1, a point also 
made by Yeh & Matzner" ('2012') . Thus one cannot simul- 



taneously have arbitrarily large arbitrarily small Pip, 
and (7^1. This issue is discussed further in Paper I. 

This problem does not affect Model GPSW, since in 
this model one ignores radiation pressure. These mod- 
els thus represent wind-dominated Hll regions, a nd are 
qualitatively similar to the models of Dopita et al. (2000) 
and [Kewley et al.| ( |2001| ). In Paper I, we show a compar - 
ison of these models with those of Dopita et al. (2000), 
and find a good match with their results. 



3.2. Dynamical effects of radiation pressure 

Having understood the effects of radiation pressure and 
winds on the internal structure of Hll regions, we are 
now ready to study their dynamical effects. 

3.2.1. Distribution of H II region radii 

In the expansion of an H II region, the radiation pres- 
sure term contributes as an additional push towards a 
faster radial expansion. To study this effect we examine 
the distribution of H II region radii produced by our pop- 
ulation synthesis code, and how it is influenced by radia- 
tion pressure. Figure [5] shows a scatter plot of the radius 
of the ionization front (Pif) versus the ionizing photon 
luminosity {S) for all the Hll regions present at one time 
step in two of our simulations, one with Pamb/^B = 10"^ 
K cm~^ (left column) and one with Pamb/^B = 10^ K 
cm"'^ (right column). We show three cases I ^trap — is 
a model where radiation pressure does not affect the dy- 
namics at all, /trap = 2 is our fiducial case, and /trap = 50 
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Fig. 5. — Radius of the ionization front Rip versus ionizing pho- 
ton luminosity S for ah the H II regions present at one time step in 
simulated galaxy (dots). We show simulations with two different 
values of Pamb/^s (10^ K cm~^ in the left column, and 10^ K 
cm~^ in the right column), and three different values of /trap (0 in 
the top row, 2 in the middle row, 50 in the bottom row). Black lines 
show the location of the stall radii for the simulations, with dashed 
lines corresponding to stalling when the pressure is gas-dominated, 
dotted lines to stalling when the pressure is radiation-dominated 
and full lines when to stalling when both radiation and gas pres- 
sure are relevant. The SFR is 1 Mq yT~^ in all the simulations 
shown, so the number of HIT regions present in each panel is ap- 
proximately the same. 



is a model where the radiation pressure is assuraed to be 
strongly trapped within the Hll region, and affects the 
dynamics much more strongly. The case /trap = cor- 
responds to H II regions that follow the classical [Spitzer 
( [l978) solution, /trap = 2 corresponds roughly to the 
va lue favored by the radiatio n-hydrodynamic simulations 
of Krumholz & Thompson pQ12a|b| ), while /trap = 50 
corresponds to t he peakof the value s adopted in the 
subgrid models of [Hopkins et al.| ( 2Qll| , where radiation 
is assumed to build up inside Mil regions and produce 
large forces. In each panel we also show with full lines 
the stalling radii, defined as the radii where the inter- 
nal pressure of the Hll region drops to Pamb- Each Hll 
region, when is created, is assigned a value of S and 
has RiF = 0. As time passes, the Hll region evolves 
and moves horizontally in the Rip versus S plane till 
it reaches this limiting line at the stall radius. Since 
Prad/Pgas dccrcascs as RiF grows at fixed depend- 
ing on the value of S and Pamb, this can happen when 



have RiF oc 5*^/^; if stalling occurs when an Hll region 
is dominated by gas pressure, then Pamb — Pgas, and 
from equation ^ we have Ri^ cx S^^^. Figure [s] shows 
also these two dependencies as dotted and dashed lines 
respectively. 

The Figure shows that radiation pressure has two dis- 
tinct effects on the dynamics. First, Hll regions with 
radiation pressure expand faster than classical ones, so 
that models are shifted to increasingly large values of 
i?iF as /trap iucrcascs. The shift from /trap = to 2 is 
relatively modest, while the gap between /trap = 2 and 
50 is somewhat larger, corresponding to nearly half a dex 
in radius. The second effect of radiation pressure is to 
increase the stalling radius. When the ambient pressure 
is small, this has a relatively small effect, because the 
stalling radius is large and most H II regions turn off be- 
fore reaching it. One the other hand, when the pressure 
is high, the stalling radius is smaller and most Hll re- 
gions stall before their driving stars evolve off the main 
sequence. In this case most H II regions are clustered up 
against the stalling radius, and the increase in stalling 
radius with /trap has very significant effects. 

3.2.2. Distribution of HII regions in the BPT diagram 

We are now ready to use our population synthesis code 
to determine where simulated galaxies lie in the BPT di- 
agram. We run three classes of models. The first, which 
we consider the most physically realistic given the ob- 
served properties of HlI regions in the local Universe, 
uses Model RPWW for the internal structures of H II re- 
gions, and uses /trap = 2 to determine their dynamical 
evolution. The other two models use /trap = (i.e. as- 
sume that Hll regions expand as classical Spitzer HlI 
regions), and use Models GPWW and GPSW for the 
internal structures. The latter choice is not fully consis- 
tent, in that with strong wind models we should use a 
wind- dominate d dynamical solution such as that of |Cas- 



Pgas < Prad, whcu Pgas > Prad, Or whcu P^g - Prad- If 

the Hll region stalls when the gas is dominated by ra- 
diation pressure, Pamb — Prad, and from equation (l9|) we 



'tor et al.'fl975"). We do not do so, however, both because 
Dopita et al. (2006b) have already explored this case, and 
"^cause observations now strongly disfavor it. 

Figure |6] shows the comparison of the three models on 
the BPT diagram for our fiducial parameter choices (see 
Table [2|. Each model represents the line ratio produced 
by summing the line emission over all the Hll regions 
present in a simulated galaxy at a given snapshot in time, 
and for each model we show 200 such snapshots, sepa- 
rated by intervals of 1 Myr. The region shown in the 
plot has been rasterized into pixels of (0.05 dex)^. The 
color in each pixel corresponds to the number of models 
that fall into that pixel, normalized by the pixel con- 
taining the most models. The plot shows several inter- 
esting results. Model GPSW, in which HlI regions' in- 
ternal structures are wind-dominated, are systematically 
shifted to lower [O lll]/H/3 and higher [N ll]/ha than the 
weak wind models. Model RPWW spans a wide range of 
par ameter space, i nclud ing some snapshots that exceed 
the 'Kewley et al. (12001) theoretical limit. These snap- 
shots tend to be immediately after the formation of a 
very large, bright, association. Model GPWW i s covers 
a smaller ra nge in the plot, and stays below the |Kewley| 



et al.| ( |2001| ) limit. 

We varied a number of the fiducial parameters, and 
found them to have little effect on the results. Param- 
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Fig. 6. — Synthetic galaxies on the BPT diagram generated with our fiducial parameters and HII region modeling. Each model is a 
single time snapshot from our simulations, with the snapshots taken at intervals of 1 Myr. The region shown has been rasterized into pixels 
of (0.05 dex)^. The color in each pixel corresponds to the number of models that fall into that pixel, normalized by the pixel containing the 
most models. The three cases shown are Model RPWW with /trap = 2, corresponding to H II regions with weak winds whose dynamics 
and internal structure include radiation pressure; Model GPWW with /trap = 0, corresponding to HII regions that are classical Stromgren 
spheres with neither radiation nor wind pressure affecting their internal structure or dynamics; and Model GPSW with /trap = 0, for which 
HII regions are wind-dominated bubbles without radiation pressure. 
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Fig. 7. — Comparison of BPT diagrams for two runs with 
-famb/^s = 10^ K cm~^ and /trap = 2 and 50. Each plus sign rep- 
resents a single time snapshot from our simulations, with the snap- 
shots taken at intervals of 1 Myr. Both runs use Model RPWW, 
and are otherwise identical to the runs shown in Figure [g] 

eters whose influence is negligible include Ma,min and 
^a,max, minimum and maximum value of the associa- 
tion mass, Ai, the Mach number used to set the width 
of the density distribution, and kp the powerlaw index 
that describes the density distribution into which Hll 
regions expand. Perhaps surprisingly, the value of /trap 
also has relatively little effect if we hold the internal mod- 
els fixed, as illustrated in Figure [71 In other words, if we 
use Model RPWW to describe the internal structure of 
Hll regions, the differences in the distributions of Hll 
region radii visible in Figure [5] as we vary /trap from to 
50 do not produce corresponaing differences in the loca- 
tions of the resulting galaxies in the BPT diagram - or at 
least the differences they produce are mostly within the 
scatter produced simply by stochastic drawing of associ- 
ation masses and surrounding densities. In Figure [l] we 
show the grid of models covering over 5 orders of mag- 



nitude both in RiF and S. The grid dramatically shrinks 
in the BPT diagram (Figure |3| causing; the small effect 
of /trap m Figure [7| Thus there does not appear to be an 
obvious way to use line ratio observations of integrated 
galaxies to measure the value of the dynamical parameter 
/trap- We stress, however, that /trap includes the influ- 
ence of wind pressure, and a wind-pressure dominated 
state can be identified, on the basis of line ratio observa- 
tions, through its effect on the internal structures of H II 
regions. In particular, wind-dominated regions cannot 
access high values of the ionization parameter and are 
limited to the lower rig ht of t he BPT diagram; see Paper 



I and Yeh & Matzner (2012) for a thorough discussion. 



3.3. Influence of galactic parameters 

While there are a number of parameters in our model 
that make very little difference to the results, the two pa- 
rameters Pamb and SFR that we use to characterize our 
galaxies do have a measurable influence. Figure |8] shows 
how the ambient pressure influences the position of sim- 
ulated galaxies on the BPT diagram. We show our three 
models computed with Pamb/^B = 10^ and 10^ K cm~^ 
(top and bottom rows). At low ambient pressure, we 
find a significantly larger spread in the simulated galax- 
ies. This is because for low ambient pressure the stalling 
radius is large, many HlI regions do not live to reach it, 
and thus Hll regions span a large range of radii. Ex- 
actly where H II regions fall in the plane of S and is 
therefore subject to a great deal of stochastic variation. 
In contrast, as show in Figure [Sj increasing the ambient 
pressure causes all the Hll regions in a galaxy to clus- 
ter along the stall radius line. In Figure [8] we can also 
see that the ambient pressure controls the overall loca- 
tion in the BPT plot, moving all the synthetic galaxies 
to a higher position in the BPT diagram, and at higher 
ionization parameter. 

Figure [9] shows the effects of varying the SFR on the 
location oi our synthetic galaxies on the BPT diagram. 
As the Figure shows, a smaller SFR leads a bigger spread 
of points in the BPT diagram. This is due to the stochas- 
tic nature of star formation at low SFRs, something that 
can also lead to large variations in absolute line fiuxes as 
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well as line ratios (iFumagalli et al. 2011 da Silva et al._ 
2012 Weisz et al.||2012[ ). If we draw a large mass for the 



next association to be created, a long time passes until 
it appears, especially when the SFR is low. During this 
phase there are no young, bright Hll regions present, 
and so the galaxy is located in the bottom-right part 
of the BPT plot. When the association finally forms, 
the galaxy's line emission becomes dominated by the re- 
sulting bright, young Hll region, which drives it to the 
top- left part of the BPT diagram. As a result, there is 
a great deal of variation in the galaxy's location. When 
the SFR is high, on the other hand, Hll regions form 
continuously, causing the population of Hll regions to 
be more numerous and uniform. We do caution that our 
mechanism for handling H II region creation may overes- 
timate the amount of stochasticity found in real galaxies, 
but that the general sense of the effect will be the same 
as we have found, even if its magnitude is overestimated. 
A more realistic formalism for handling the problem of 
drawing association masses and birth times subject to 
an overall constraint on the star formation rate is imple- 
mented in the SLUG code ( da Silva et al.||2QT2 ); adding 
this formalism to our code is left for future work. 

4. COMPARISON TO OBSERVATIONS 

Having understood the physics that drives the location 
of galaxies in the BPT diagram, we are now in a position 
to compare our models to observations. Such observa- 
tions come in two varieties: spatially resolved ones of 
individual Hll regions or portions of galaxies, and un- 
resolved ones in which the line fluxes from all the Hll 
regions in a galaxy are summed. Since our code produces 
collections of stochastically-sampled Hll regions, we can 
compare to both. For reference and to facilitate compar- 
ison, we show in both cases unresolved observations of 



local galaxies from th e SDSS ( |Brinchmann et al. 2004 
Tremonti et al.|["2004| alo ng with a tit to this sequence 



( [Brinchmann et al. J2008 |, the empirically-determined 
line separati ng star-forming galaxies from AGN (|Kauff-| 
|mann et al.|[2003), and the Kewley et al.| (|200ip theo- 
retical upper limit to star forming galaxies. The single 
Hll region sequence and the SDSS star forming galaxy 
sequence overlap, at least in the upper left part of the 
BPT diagram, while high redshift galaxies seem to ere 
ate a different sequence, upward and to the right ( Liu 



et al. "20081 'B rinchmann et al.|[2QQ8l [Hainline etHJIm^ 
Erb et al. 2010^ 

Figure [lO] shows a collection of observations of single 
Hll regions and pixel by pixel observations taken from 
the literature. For the local Universe, our comparison 
data s et consists of single Hll region s from NGC 136 5 
( |R8^fc Wals h 1997), NGC 1313 (W alsh fcR2 ^ p97|) 
and t he Orion region in our own galaxy ( [Sanchez et al. 
|2007l). We also plot individual pixels in three lensed 
galaxy at z ~ 2 from [Jones et al. (2012), which scatter 
about a locus that passes close to the loc ation of Orion 
in the BPT diagram. As pointed out by [Walter et al. 
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Fig. 10. — A comparison between resolved HII regions and in- 
dividual H II regions produced in our simulations in the BPT dia- 
gram. We show H II regions in galaxies at z = with blue symbols; 
the galaxies shown are NGC 1365 (Roy &; Walsh 1997) (triangles), 
NGC 1313 (Walsh & Roy 1997) (diamonds), and the integrated 
value for the Urion nebula ( ISanchez et al.,^2007| (star). We also 
show single spatial pixels measured wit h UdIHIS for thre e lensed 
galaxies at redshift z ~ 2.0 — 2.4 from ([Jones et al.[[2012|) ( green, 
red, and orange). Contours represent galaxies trom the bDSS cat- 
alogue (Brinchmann et al. 2004), enclosing respectively, 5, 10, 20, 
50, 90 and 99% of all galaxies in which the four emission lines are 
detected at a greater than 3a significance level. Blue and purple 
points show the results of our models using fiducial parameters and 
Pamh/kB = 10^ and 10^ K cm""^, respectively. Finally, the dot- 
ted line is the theoretical upper limit of Kewley et al. (2001), the 
dashed line is the empirical AGN - star-forming galaxy separating 
line (^Kauffmann et al.^2003J, and the dot-das hed line is the fit to 
the SUbb sample ot [ljrinclLmann et al.[ ( |2008t . 

regions produced in two different simulations, one with 
Pamh/ks = 10^ K cm~^ one with 10^ K cm~^. These two 
cases should roughly bracket what we expect for Milky 
Way-like galaxies and for the dense, more strongly star- 
forming galaxies found at high redshift. The plot shows 
that our simulations are able to roughly reproduce the 
locus of observed H II regions in the BPT diagram for a 
reasonable range of ambient pressures. We cannot repro- 
duce most of the Hll regions in NGC 1313, because the 
galaxy has a metallicity lower than solar and our model 
considers only solar metallicities. Lower metallicity pro- 
duces a shift of the m odels towards lower [N ll]/Ha values 
( Dopita et al.[[20QQ ). The pixel by pixel high-z galaxies 
are best tit by the models with high Pamb, consistent 
with observations that these galaxies have high surface 
and volume densities. 



Figure pT] shows the comparison with integrated galaxy 
measurements; these come from the SDSS for the local 
Universe, and from a variety of surveys at high- 2:. Many 
SDSS star forming galaxies lie in the lower part of the 
star forming sequence due to the presence of a diffuse 
warm comp o nent in the interstellar medium. jBrinch- 



(^2009), the SFR surface density of Orion is similar to "mann et alj ( 2004 ) point out that a significant amount 



that of a high redshift object undergoing a burst of star 
formation. 

On top of these data, we overlay the results of our sim- 
ulations using model RPWW, which we consider most 
realistic based on observations of nearby Hll regions. 
The results shown are single snapshots of all the Hll 



of the emission line flux in these galaxies comes from 
the diffuse ionized gas, rather than from Hll regions. 
The combination of the diffuse ionized gas and the Hll 
regions typically has a lower effective ionization param- 
eter, and compared to Hll regions alone it shows an 
enhanced [N ll]/Ha and depressed [O lll]/H/3 ( ,Mathis[ 
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2000). Therefore, we only expect our models, which do 



not include the diffuse ionized gas, to reproduce the up- 
per part of the star forming sequence of the SDSS. 

In Figure [TT] we also overplot the whole-galaxy results 
produced by our code. As the plot shows, while we are 
able to reproduce the full spread of individual Hll re- 
gions, our simulations of whole galaxies cover a more 
limited range of BPT than the observations. In partic- 
ular, we tend to underpredict the observed [N ll]/Ha 
ratios. There are several possible explanations for why 
we might successfully reproduce individual Hll regions, 
even in high-z galaxies, but not fully cover the range of 
integrated galaxy properties. One we have already dis- 
cussed in the introduction: the offset at high-z may be 
due to the contribution of a weak AGN, which our mod- 
els obviously do not include. A second possibility is that 
the contribution of diffuse ionized gas to the line ratios 
cannot be neglected even in these high redshift galaxies. 
Another is that our weighting of the different HlI re- 
gions is incorrect because the association mass function 
is different than the —2 powerlaw we have adopted based 
on local observations, or because of biases introduced by 
dust extinction, despite the extin ction-independ ent na- 
ture of the BPT line ratios (see |Yeh fc Matzner||2012| . 
A fourth possibility is that our lognormal distribution of 
densities provides a poor fit to the true range of densi- 
ties into which H II regions expand in high-z galaxies, so 
that the amount of time individual H II regions spend in 
the upper left versus the lower right parts of the BPT 
diagram is off in our models. 

As a last possibility, we recognize that the ability of 
our models for individual HlI regions to reproduce the 
observations of Orion very well (Figure 10) may be partly 
a matter of good luck. Our models are not designed to 
mimic the champagne flow phase of young (< 10^ yr), 
compact Hll regions. In particular, we assume a state 
of quasi-static force balance wh ich holds only approxi- 
mately in accelerating flows; see Yeh & Matzner (12012) 



3.4 on this point. Indeed, Orion does not resemble the 
typical HlI region - e.g., a few million years old and at 
the stalling radius - in our galaxy simulations. It is pos- 
sible that the different distribution of high-z galaxies in 
the BPT plot as compared to local SDSS galaxies is due 
to the higher pressure environment in the former, which 
keeps the Hll regions longer in a champagne flow-like 
phase. Future studies might assess our models' accuracy 
in the champagne phase, extend their range of validity, 
and quantify the importance of this dynamical detail for 
high- 2; galaxies. 

5. DISCUSSION AND CONCLUSIONS 

Motivated by recent observations suggesting that H II 
regions are shaped much less than expected by the pres- 
sure of shocked stellar wind ^as , and much more by direct 
radiation pressure ([Ha rper-Cl ark fc Murray||2009 Lopez 
et al.|2011 Yeh fc ]CTa!;zner|20l^ ) , we revisit the problem 
ot determining the line tlux emitted by a population of 
Hll regions. We adopt as our default a model of Hll 
regions where the pressure of winds is subdominant, and 
radiation pressure is not neglected, and we compare this 
result to traditional models with strong winds and weak 
radiation pressure. In Paper I we discuss how we gener- 
ate grids of static, single Hll regions, with a wide range 
of sizes and ionizing luminosities, with varying strengths 
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Fig. 11. — A comparison between observed unresolved galaxies 
and simulated produced by our code in the BP T diagram. Con- 
tours and lines are the same as in Figure |10| Diamonds show 
galaxies at z ~ 1 — 1.4 taken from the DEIEP2 survey Shap-J 
[ley et a l. (2005); |Liu et a l. (2008); triangles show galaxies at 
z ^ 2 taken from IhJrb et al. (2UU6 2010) ; lensed galaxies at a 
varie tv of re dshifts and samples are indicat ed bv inverted trian - 
gles dLe moine -Busserolle et al. 2003|), circ les ([Hainline et al.|20Q9| , 
squares ( [Kicliard et al. 2U11), and stars ( Jones et al. 20 vA. hJlue 
and purple plus signs show the results for integrated galaxies in 
of our models using fiducial parameters, Model RPWW, and am- 
bient pressures of Pamh/ks = lO'^ and 10^ K cm-^. Each point 
represents a single time snapshot taken at 1 Myr intervals from a 
simulation that evolves for 200 Myr. 

of winds and radiation pressure. In this paper we con- 
struct dynamical expansion models for these H II regions, 
and explore how changing the strength of winds and ra- 
diation pressure affects their line ratios in the BPT dia- 
gram. We find that radiation pressure has two important 
effects. First, Prad changes the internal structure of the 
H II region, creating a density gradient towards the outer 
shell. This affects the expected line emission, allowing 
the H II region s to exceed the upper limit form starburst 
models set by ^Kewley et al. (2001). Second, radiation 
pressure provides an extra boost to the expansion, lead- 
ing to larger radii at earlier times. 

We embed these models in a population synthesis 
code that generates galactic collections of stochastically- 
generated Hll regions expanding into a turbulent 
medium. The code follows H II regions as they are born, 
evolve, stall and die. Using this code we predict the in- 
tegrated line emission of galaxies as a function of several 
galactic properties. We find that the two most important 
ones in controlling where galaxies appear in the BPT di- 
agram are the ambient pressure, which shifts galaxies up 
and to the left as it increases, and the star formation 
rate, which affects the amount of stochastic scatter in a 
galaxy's line ratios. 

We compare with observations in two distinct ways. 
First, we select single Hll regions observed in the local 
universe and pixel by pixel observations of z ^ 2 galaxies, 
and we compare these to the distributions of individual 
Hll regions produced in our model. We show that our 
model produces good agreement with the observations 
for reasonable ranges of SFR and ambient pressure. The 
high redshift pixel data are best reproduced by Hll re- 
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gions evolving in a high pressure medium and with high 
SFR, which we interpret as a sign of intense star forma- 
tion in a dense interstehar medium, consiste nt with the 



observed properties of high-z galaxies (e.g. [Elmegreen 
|et al. 2009j [Genzel et aL] |2Qll). 



Second, we compare mtegrated galaxies from the SDSS 
catalogue and the high redshift universe to our synthetic 
galaxies. We find that, while we are able to reproduce 
the spread of individual Hll regions, our models for in- 
tegrated galaxies cluster too tightly compared to the ob- 
served range of line ratios in real galaxies, particularly 
at high-z. This might be due to a number of factors. 
One possibility is that the lognormal distribution of the 
ambient density we have adopted is a poor description of 
the density distribution in high-z galactic disks. Another 
possibility is that winds might be important at high red- 
shift or that the presence of the diffuse ionized medium 
is not negligible. A third possibility is that a higher pres- 
sure environment in high-z galaxies keeps the H II regions 



longer in a champagne fiow-like phase. One last possi- 
bility is that z 2 star forming galaxies may contain an 
AGN that partially contributes to the line emission. We 
leave these possibilities as a subject for future work. 
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